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We develop a systematic expansion method of physical quantities for the SK model and 
the finite-dimensional ±J model of spin glasses in non-equilibrium states. The dynamical 
probability distribution function is derived from the master equation using a high temperature 
expansion. We calculate the expectation values of physical quantities from the dynamical 
probability distribution function. The theoretical curves show satisfactory agreement with 
Monte Carlo simulation results in the appropriate temperature and time regions. A comparison 
is made with the results of a dynamics theory by Coolen, Laughton and Sherrington. 
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§1. Introduction 

Dynamics plays essential roles in understanding experimental observations on spin glasses. It is 
very difficult, however, to develop systematic theoretical methods to investigate dynamical behavior 
of spin glass models, and numerical approaches have been the major source of information on 
dynamical properties of finite-ranged spin glasses until recently. Recent activities in analytical 
studies on the dynamics of the Sherrington-Kirkpatrick (SK)0) model of spin glasses include closed- 
form solutions for dynamical correlation functionsll' and construction of evolution equations for 
single-time macroscopic quantities .SS 

In the latter approach, Coolen, Laughton and Sherrington (CLS) derived a closed-form evolution 
equation for the averaged single-site spin-field distribution function under a few assumptions on 
microscopic properties of the system. The averaged spin-field distribution function thus obtained 
provides sufficient information to determine the time evolution of single-time macroscopic quan- 
tities. It gives quite good agreement with the evolution of macroscopic quantities obtained by 
computer simulation of the microdynamics, although a high-temperature expansion of the dynam- 
ical microscopic probability distribution!) involved quantities not expressible solely in terms of the 
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spin-field distribution of CLS. 

Tlie purpose of the present paper is to further develop the high-temperature expansion technique!^ 
and derive explicit expressions of physical quantities to the third order in the inverse temperature. 
The obtained series results are still too short to discuss critical properties around the spin glass 
transition temperature. However, our method has an advantage that it works not only for the 
infinite-range SK model but also for finite-dimensional systems, in particular the nearest neighbor 
lb J model. Thus the present approach is a first step toward a systematic investigation of dynamics 
of finite-dimensional spin glasses. 

In the next section the problem and the method are formulated. The explicit expression of 
the dynamical probability distribution function is derived there to the third order of the inverse 
temperature. A method to evaluate expectation values of physical quantities using the dynamical 
probability distribution function is formulated in §3. This formulation is used in §4 to determine 
the coefficients of series expansions of several physical quantities. In §5 the averaged spin-field 
distribution function is evaluated by the expansion method. The results in these two sections are 
used in §6 to discuss the limit of applicability of the CLS theory. The final section is devoted to 
general discussions. 

§2. Dynamical probability distribution function 

The essence of our theory is a high-temperature expansion of the microscopic probability dis- 
tribution function. The basic formulation of this expansion was developed in Ref. |5| for the SK 
model of spin glasses. We describe here briefly the method generalized to include finite-dimensional 
models. 

The model we consider is an Ising spin glass with Hamiltonian 



where <t denotes the spin configuration and the {Jij} are quenched randomly distributed; {ij) 
denotes a pair of sites. The probability distribution function pt{cr) obeys the master equation: 



where is a single spin fiip operator, F^^^cr) = ^{ai, • • • , —a^, • • • , ctn), being the system size, 
and we use a transition rate of the heat bath method 





{ij) 




(2.2) 




(2.3) 



Here /? is the inverse temperature and hk{cr) denotes the local field on site k. 
Let us solve the master equation (|2.2D by a high-temperature expansion in the form 




(2.4) 
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Inserting eq. (2.4) into the master equation ( p.2| ) and expanding the result in powers of P, we obtain 
differential equations for the functions appearing in the exponent of eq. (|2.4|): 



^ = ^Y.i^kft + 2akhk) (2.5) 

k 

^ = ^ E {li^kftf + Akgt + cJkhkAkf^j (2.6) 
^ = ^ E I ^(^fc/*)' + ^kftA^gt + AkUt 

+ Gkhk ( Afc/t)2 + Afc j - ^fjfc/i^ I , (2.7) 

where A^/t = ft{Fkcr) - ft{a). 

We now restrict ourselves to the nearest neighbor it J model on a finite-dimensional lattice. The 
SK model is recovered by taking an appropriate limit of infinite coordination number. The system 
is assumed to be initially in equilibrium at the inverse temperature Pq: 

Pt=o{cr)=exp[-poH{cT)]. (2.8) 

Equation (2.5) is easily solved to yield 

McT)=a{t)H{cT), (2.9) 

where 

a{t) = ae~2* - 1. (2.10) 

The parameter a denotes 1 — Pq/ (3. It is necessary that (3q and (3 are smaller than (3c, the inverse 
of the critical temperature, because the high-temperature expansion is valid only when the system 
stays in the paramagnetic phase. 
Using eq. ( |2.9| ), the next order differential equation (2.6) is written explicitly as 

^ = ^ E ^kgt + (a^e-^* - ae-2*) ^ hi (2.11) 

k k 

This equation has the following solution 

gt{cT) = h{t)Y,hj + b2{t), (2.12) 

i 

where 

„2 / „,2 



h{t) = -^e-'^-{at-^\e-^^ (2.13) 



and 



b,it) = NzJ' { -e-^* + at - - + - e-^' + T " 2 ^ ^^"^^^ 
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with z being the coordination number. Since eq. (2.14) is independent on the spin configuration, 
this term does not affect the fohowing argument and shah be ignored hereafter. 



Similarly, using eqs. (|2.9| ), (|2.10| ), ( p.l2|) and ( 2.13 ) for and gt, eq. (2.7) is rewritten as 

AUi 1 



dt 



AkUt - 4zj2(2ae-2* - l)bi{t)H{a) 



k 



-2{2ae-^' -l)bi{t)Y,Jkihkhi. 

k,l 

The solution of this equation is 

ut{cT) = ci{t)H{a) + C2{t)Yaihi + cs^t) ^ Jijhihj, 

i 

where the coefficients satisfy 

ci{t) = -2ci(t) - (12z - 8) J2c2(t) 
-4zj26i(t)(2ae~2* -1) 



(2.15) 



(2.16) 



C2{t) = -4c2(t) + -a^e - a~e 
C3(t) = -2c3(t)-26i(t)(2ae 



(2.17) 



-2t 



1). 



In deriving these equations we have assumed that there are no triangular loops on the lattice in 
the sense 

JklJlmJ-mk ~ 0; 

where indices denote neighboring sites. This assumption excludes several types of lattices from 
our considerations such as the triangular lattice. Although it is possible to remove this simplifying 
assumption, resulting complicated formulas do not lead to new physics. 
Equations (2.17) are solved under the initial condition ci = C2 = C3 = as 



ci(t) = J' 



[Az - 2)a^ - (lOz - 4)Q^t [.e~^* + <! 2zaH 



2zar 



2z-l)a3 + (4z-2)a2 e-2t 



C2(t) = -3«'(e-'* 



1 



-4i 



C3(t) = --a3e-6* + {a'-- 



a^-2aH]e-^' 



-a +-a)e . 



(2.18) 
(2.19) 

(2.20) 
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Since eqs. (2.5), (2.6) and (2.7) are first-order differential equations, the above expressions ( p.9D , 
( p. 12 ) and (2.16) with coefficients obtained above represent the unique solution. We have therefore 
obtained the dynamical probability distribution function to the third-order in [3 for the it J model 
as 



Pt{c7)= exp 



l3a{t)H{a)+P^hi{t)Y,hj 

i 

+ (3^(^ci{t)H{a) + C2{t)J2^ihf 
+ csit) X] Jijhihj j H 

i,j ) 



(2.21) 



where a(i), 6i(i), ci(t), C2(i) and C3(t) are given in eqs. (^l0|) , (|2l3|) and (2.18) - (2.20). 

We have considered the it J model on a finite-dimensional lattice with z nearest neighbors. The 
results derived previously for the SK modefi are recovered by taking the limit N ^ oo with 
z = N-landJ = J/VN (J ~ 0(1)). 



The above solution (2.21) is invariant under the gauge transformation aiTi,Ji. 



[Ti = ±1). This invariance is naturally expected from gauge invariance of the master equation 
(p.2|). In general, any term in the expansion can be expressed in terms of lij = Jijaiaj which is the 
basic building block of all gauge invariant quantities. More explicitly, terms to the third order are 
written as 



{ij) 



ilm 



^ ^ ^ihi — ^ ^ lillimlin 

i ilmn 

^ ^ Jijhihj — ^ ^ lijiiiljmi 
ijlm 



(2.22) 

(2.23) 
(2.24) 
(2.25) 



where all the indices of the Fs stand for nearest neighbor sites. Higher-order terms should be able 
to be written similarly. 

The series-expansion solution obtained in the present section forms the basis of estimation of 
physical quantities in the following sections. 

§3. High-temperature expansion of physical quantities 

We now apply the series expansion solution for the probability distribution function to the ex- 
pectation value of a physical quantity O: 

J2cT0{<^,{Jij})Pt{(^) 



0{a,{J,j}))^^ 



EcrPt(o-) 



(3.1) 
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The first-order solution (|2.9| ) with eq. ( p.lO|) gives 

pt{cr) = exp - ae''^*)H{a) = exp[-/?efri^(cr) 



(3.2) 



This equation shows that the dynamical expectation value (3.1) can be calculated using a time- 
dependent effective inverse temperature /^cfr- Therefore we may regard the system as being in 
equilibrium at temperature T/{1 — ae~^*) at any given time t. The expectation value with respect 



to this effective Boltzmann factor (3.2) will be denoted by (• • ■)i in the following. 
The higher-order terms in the series expansion are taken into account by the expansion 

exp|/?a(t)/7 + pHiit) hi + PHci{t)H + • • oj ^ 

The thermodynamic average of a physical quantity 0(cr, {Jij}) then becomes 
, , (0 + pHi{t)0j:ih^ + PHci{t)0H{(T) + 

(0{(T,{Jij})) ) 

^ {l + f3%l{t)E^hf + P^{c,{t)H{a) + •• 

- (O), + /3%(t)((0E^?>, - (0)XE^?>,) +/3^|ci(t)((0if(-))^ - {0)^{H{a))^^ 



(3.3) 



+ 



(3.4) 



The configurational average of this equation is expressed as 



O 



O 



+ 



(3.5) 



i 

+ f\cx it) Gov (0H{(t)) ^ + C2 it) Gov (o ^ ^^i^^) ^ 

+ C3(t)Gov(0^ Jij/li/lj)^ I. 



Here (0)„ stands for the thermodynamic average obtained by truncating the exponent of eq. ( |2.4 
at the nth order in /?, and Cov(^i3)i = I^AB)\ — {A)i{B)i represents the covariance of A and 
B at the effective inverse temperature PeS- 

Generally it is quite difficult to evaluate the expectation value (0)i and covariances on the right- 
hand side of eq. (3.5) for the finite-dimensional itJ model. An exception is the first-order term 
(O) 1 evaluated for a gauge-invariant quantity O on the Nishimori line defined by tanh PJ = 2p—l, 
where p is the probability for Jij to be positive.!^ In the case of the SK model, we can derive explicit 
expressions for all terms in eq. (3.5) by the replica method. 
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Therefore we have to mtroduce an approximation to proceed further for the it J model. For this 
purpose, we first focus our attention on the equilibrium expectation values of terms appearing in 
the expansion (3.5) evaluated on the Nishimori line: 



H(a] 



i 



1 



iVJztanh/3J 



iVj2 + z(z - l)(tanh/3J)2} 
= NJ^{z{?,z - 2)tanh/3J 



(3.6) 
(3.7) 



+ 



z(z-l)(2-2)(tanh/3J)3} 



+ z(z-l] 



NJ^{z{2z-l)tanh(3J 



'(tanh/3J)3} , 



(3.8) 



(3.9) 



where (•••)e denotes the equilibrium expectation value. The first of these equations, eq. (3.6), 
indicates that the factor tanh /? J can be replaced by the expectation value of the Hamiltonian. For 
instance. 



NJ^z - 2J{z - 1) tanh pj (H{(t) 



(3.10) 



In equilibrium, both sides of eq. ( 3.10| ) are of the order of N, and fluctuations around the expectation 
value are of order \/iV- Thus, given a typical equilibrium spin configuration, the value of the 
quantity 

i 

is almost certainly equal to that of 

Nj'^z - 2J{z - 1) {tanh (3 J) H (a). 

This argument does not apply to temperatures near the critical point where fluctuations are not 
negligible. 

The situation should not change drastically if the system is close to equilibrium if not in equilib- 
rium precisely. Thus the following replacement may be a reasonable approximation if the system 
is not far from equilibrium. 



N.fz - 2J{z - l)(tanh/3J)i?(cr). 



Further a posteriori justification of this approximate replacement comes from agreement of the 
resulting values with simulations as shown in the following. 
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We therefore apply similar replacements to all terms in the exponent of eq. (2.21) using the 
relations (3.7) to (3.9). The result is 

Pt{(T) ~ exp(^Pa{t)H{a) + ^^bi{t)i^NJ^z - 2J{z - l)(tanh/3J)i7(cr)| 
+ /3^|ci(f) - 2c2(t)J^ (^3z-2 + {z- l){z - 2)(tanh/3J)2j 

- 2c3(t) {2z-l + {z- l)2(tanh pjf^ |i^(cr)^ . (3.11) 

This approximation is valid for a system near equilibrium on the Nishimori line. 

This result (3.11) indicates that the dynamical probability distribution is approximated by the 
equilibrium Boltzmann factor with effective inverse temperature 

(/3eJ)±J = -Mt) + 2(3%i{t)J{z - l)tanh;gj - /J^ x 

{ci(t) - 2c2(t)J^(3z - 2) - 2c3(t)J^(2z - 1)} . (3.12) 

In the case of the SK model, the thermal expectation values on the left-hand side of eqs. (3.6) - 
(3.9) can be evaluated using the replica method. The corresponding (3^^ for the SK model is then 
obtained as 

(/3S)sK = -Pa{t) + 13^ {2phi{t) - ci{t) 
+ 6J^C2{t)+4:J^cs{t)), (3.13) 

where J'^/N represents the variance of distribution of Jij. The above formula for the SK model can 
also be derived from eq. (3.12), which is valid on the Nishimori line, by taking the limit N — > oo, 
z = N — 1 and J = J/\//V and applying the approximation tanh/?J ~ pj in the it J model. Note 
that equilibrium properties of the SK model in the paramagnetic phase are independent of Jq, the 
center of distribution of Jjj.B Therefore the results valid on the Nishimori line, eqs. (3.6) - (3.9), 
which passes through the paramagnetic phase toward the multicritical point, remain true in the 
whole region of the paramagnetic phase of the SK model. This is the reason why eq. (3.13) is valid 
in the whole paramagnetic phase of the SK model whereas eq. (3.12) for the finite-dimensional 
nearest-neighbor it J Ising model can be used only on the Nishimori line. 

In the form exhibited in eq. (3.5) the direct expansion is useful in practice only for the infinite- 
ranged SK model or to the level of the first-order term. Quantitative estimation of validity of these 
methods will be studied in the next section. 

§4. Explicit evaluation of expansion terms 

We now apply the formulas derived in the previous section to the following physical quantities, 
H, J2i hf, J2i (^i^i and J2ij Jijhihj. These four quantities constitute the expansion of the dynamical 
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probability distribution function to the third order in P as in eq. (2.21). One of the purposes to 
choose these quantities for a test ground of the method developed in the previous section is to 
estimate the order of magnitude of various terms in eq. (2.21) for consistency check of the theory 
of Coolen, Laughton and SherringtonBS as discussed in §6. 

The first-order approximation for the it J model is obtained by replacing /? by PeS in eqs. (3.6), 
(3.7), (3.8) and (3.9). The corresponding expressions for the SK model are found from those of the 
lb J model by taking the limit of infinite dimensionality, or alternatively, by applying the replica 
method as described in Appendix A. The results are 



L ij J 

The third-order formula (3.5) can be explicitly evaluated only for the SK model as mentioned 
before. The necessary covariances are calculated in Appendix A using the replica method: 




(4.1) 



(4.2) 



(4.3) 



(4.4) 
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\ i ij / 



= iVj6 (12 + + 18 J^^) 

Inserting these results into eq. (3.5), we obtain 



1 



+N[ -J'ci(t) - 3J*C2(t) - 2j4c3(t) 



E 



E 



+ 2NJ%i{t)l3^ 

[- l8NJ^bi{t)PesP^ 



+iV-3J^ci(t) + 247^2(0 + 12rc3(t)/?' 



E Jijhih^ 



L 2J 

X/ 



+ l2NJ%i{t) 



^E '^■i-j^i^j/^ 

;/?eff/?^ + iv(-2j4ci(t) + 12j6c2(t) + 10j6c3(t))/3l 

The approximation of a Boltzmann form with the inverse temperature, eqs. (3.12) and (3.13), is 
also used in the following analysis. 

We have compared the theoretical curves obtained by the approximations (• • •)!, (• • •)3 and 
with each other and with Monte Carlo simulations for the SK model. Figures 1 to 4 show the results 
for H, J2i hf, J2i '^if^i aiid J2ij Jijhihj, respectively. The initial temperature was Tq = oo and the 
final temperatures were set to T = 2 and 5 in each figure. The system size of simulations was 
iV = 1000 with 50-sample averages and standard deviations indicated. The center of distribution 
is Jo = 0, and the temperature is expressed in units of J/Zcb or J /k^. 

For all quantities the theoretical curves follow simulations faithfully at relatively high tempera- 
ture, T = 5. When the temperature is lower, T = 2, the first-order approximation deviates in the 
intermediate time region, t = 1 ~ 2, from simulations as well as from other theoretical curves. 

We next investigate the nearest neighbor it J model. The straightforward expansion (3.5) cannot 
be used in this case because the covariances are difficult to evaluate explicitly. We therefore compare 
the first-order term and the -approximation with simulations. 

The results for four quantities of the it J model on the square lattice are shown in Figs. 5 to 8. 
We have set Tq = oo and the final temperatures were T = 2 and 5. Monte Carlo simulations were 
carried out for the system size = 50^ on the Nishimori line. The plots show results averaged 
over 50 samples with standard deviations indicated. 

The critical temperature is Tc ~ 0.96 for the two-dimensional it J model on the Nishimori line,i^ 
while the SK model has Tc = l.^ In addition, the two-dimensional model has the problem of 
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Griffiths singularity at Tq = 2.270 Therefore it is not surprising that the T = 2 results do not 
reproduce simulation data very well while the T = 5 case does. The slightly better situation at 
T = 2 for the SK model than in the it J model may be related to the absence of Griffiths singularities 
in the SK model. 

§5. Spin-field distribution function 

We now calculate the dynamical single-site spin-field distribution function by high-temperature 
series expansion. This distribution function has more information than individual physical quan- 
tities in the sense that the physical quantities such as [{H{cr))], [(^ihf)] and [(^iCTih^)] can be 
derived from the distribution function. It should be noted that Jijhihj)] cannot be calculated 
from the present distribution function, because this distribution function does not carry information 
on correlations between fields at different sites. 

The averaged single-site spin-field distribution function is defined by 



Pih,a) 



(5.1) 



where cr = ±1 should not be confused with the microscopic spin configuration of the whole system 
<T. The equilibrium spin-field distribution has been calculated by Laughton et a/.i^ as 

(/i-aj2/52)\ 



P{h,a) 



1 



2V2ttJ 



- exp 



2J2 



(5.2) 



Now we consider the dynamical single-site spin-field distribution function. This function should 
reduce to eq. ( |5.2| ) in the equilibrium limit. It is convenient to introduce the characteristic function 
(Fourier transform): 



Pih,a) 



dk 
2^ 



G{k, a) exp(— ifc/i), 



or 



G{k,a) 



n' 



(5.3) 



(5.4) 



The method of high-temperature series expansion developed in the previous sections is applicable 
to the evaluation of eq. (^.4]). Detailed calculations are described in Appendix B. The result is 



Oik, a) 



1 



2 



k'J 



■ exp 



|cos(A;J^/9efr) + icr sin (A:J^/3efr) | 
+ /J^bil - AkJ'^peS sm{kJ^pes) - k^J^ cos{kJ^/3es) + i(T (4A;J^/3efr cos{kP/3e{i) - k^J^ sm{kJ^pes)) } 
+ l3^i^ci{kPsm{kJ^/3es) + i(T (^-kJ^ cos{kJ^ /3es)) }+C2{ - 3k^ cos{kJ^ /S^s) 
- {J^{'ok0l^ - k^) + QkJ^) sm{kJ^P^s) 

+ ia (-3A:2 JVeff sm{kJ^ M + {J^{6kp'^fi - k^) + 6kJ^) cos{kP(3es)) j+csj - 2^;^ J^eff cos{kJ^ M 
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+ ia (^-2k^J^Pes sm(A:P/3eff) + {QkJ^p^s + 4A: J^) cos(A:P/3eff)) } 
Substituting eq. (5.5) into eq. (|5.3|), we arrive at the expression of the spin-field distribution: 



(5. 



Pih,a) 



1 



2x/2ttJ 



+ P^ci{-aih-aJ^P,s) 



X exp 



{h - aJ^M' 
2J2 



(5.6) 



It is straightforward to check that the expansion results of [{H{cr))], [{^{hf)] and [(X^i^j^f)] 
given in §4 are recovered from eq. (5.6) by appropriate integrations. 

In equilibrium the distribution (5.6) becomes Gaussian. However it is not Gaussian in general 
even in the paramagnetic phase. As the temperature is lowered, the distribution deviates from 
Gaussian more significantly. 

The results of the analytical calculation and computer simulations are compared in Fig. 9. 
Simulations were performed for N = 5000 with 100 samples at T = 2. We may conclude that the 
third-order result (5.6) agrees with simulations in a rather satisfactory manner at this temperature. 

§6. Discussions on the CLS theory 

Coolen, Laughton and Sherrington (CLS)i'& developed a theory of dynamics of the SK model. 
They derived the closed-form evolution equation of a physical quantity which determines macro- 
scopic behavior of the system. An important assumption in their derivation was that the average 
of a macroscopic physical quantity can be calculated under the ansatz of equipartitioning of the 
dynamical probability distribution function. Here equipartitioning means that the probability dis- 
tribution function is constant once the value of the single-site spin-field distribution function P{h, a) 
is given. 

The existence of the correlation term of fields at different sites, Jijhihj, in the probability 
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distribution function (2.21) shows that the equipartitioning ansatz of CLS is not generahy true and 
suggests that the CLS theory is an approximate one-S We investigate in the present section the 
effects of the field-correlation term on physical quantities to determine the degree of applicability 
of the CLS theory. 

Figure 10 shows the third-order series-expansion curves and simulation results, superimposed for 
comparison, of various physical quantities. This figure indicates that the field-correlation term is not 
necessarily small compared with other terms appearing in the series expansion. The time evolution 
of the correlation coefficients of Jijhihj and the other quantities is shown in Fig. 11 obtained by 

(3) 

the /3gg -approximation. The correlation coefficient between A{cr, {Jij}) and Jijhihj is defined 
as 



Figure 11 shows that J2ij Jijhihj is closely correlated with H and X^i '^j^f but is not with hf in 
the initial time region. Since the absolute value of correlation coefficients with H and <^j^f is 
close to unity, one may suppose that the field-correlation term ^j^- Jijhihj is approximately replaced 
by H or J^i'^ihf, which may be taken as a support for the CLS theory as an approximation. We 
can also obtain results similar to Fig. 10 for the nearest neighbor it J model. 

It may be useful to see the behavior of various terms in the series expansion with the time- 
dependent coefficients a{t),bi{t), ■ ■ ■ , 03(1) taken into account as appearing in eq. (2.21). Third- 
order series estimates are given for the SK model in Figs. 12 and 13 and in Fig. 15 for the 
two-dimensional ibJ model in the high temperature region. Simulation results run almost along 
the same curves and are omitted in these figures. Only simulation results are shown at lower 
temperatures, T = 1 for the SK model and T = 2 for the it J model, because the expansion breaks 
down at these temperatures; see Figs. 14 and 16. 

Our observation is that the field-correlation term multiplied by the coefficient, /3^C3(t) J2ij Jijhihj, 
is very small compared to first and second order terms in the high-temperature region, T = 5, both 
for the SK and itJ models. However, this is not necessarily the case in intermediate time regions 
at lower temperatures; see Figs. 14 and 16 which show effects around (t ~ 2). Naively one might 
therefore expect to need to choose appropriate parameter regions when applying the CLS theory. 
However, trends of strong correlations for J2ij Jijhihj with H and J2i '^ihf as seen in Fig. 11 allow 
J2ij Jijhihj to be effectively approximated by H or J2i'^ihi, giving a support for the CLS theory 
as an approximation irrespective of the size of the term J2ij Jijhihj. 




(6.1) 
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It is useful to compare the macroscopic spin-field distribution P{h, a) obtained by the CLS theory 
and that by the high-temperature expansion as shown in Figs. 17 and 18. Figure 17 shows the 
first-order cumulant, or the average, of P{h,l),P{h,—l) and P{h,l) + P{h,—1) at T = 2. The 
simulation data are also displayed which do not seem to give a clear advantage either to the CLS 
theory or to the high-temperature expansion. The second-order cumulants are shown in Fig. 18. 
The upper set of curves are for P{h, 1) -|- P{h, —1) and the lower curves correspond to P{h, 1) and 
P{h, —1). This figure indicates that the simulation results follow the CLS theory more closely than 
the high-temperature expansion in this temperature region. Therefore the CLS theory serves as 
quite a good approximation for the macroscopic spin-field distribution function. 

§7. Discussions 

We have derived explicitly the time dependence of macroscopic physical quantities for the SK 
model and the nearest neighbor itJ model on the Nishimori line. Our formulation gives the 
microscopic dynamical probability distribution function in the form of a high-temperature series 
expansion. Evaluation of physical quantities has then been reduced to calculations of various 
expectation values in equilibrium. Such expectation values can be calculated to any order in the 
inverse temperature P in the case of the SK model. For the ± J model on a finite-dimensional lattice, 
the situation is more complicated and explicit evaluation is possible only to the first order in /3 on 
the Nishimori line. Wc therefore introduced an approximation in which the dynamical probability 
distribution has the same form as the equilibrium Boltzmann factor with a time-dependent effective 
temperature. 

The resulting expressions show excellent agreement with numerical simulations at high tempera- 
tures. Deviations are observed at lower temperatures in the intermediate time region. Wc anyway 
think it a significant progress that a systematic method to evaluate explicitly the time dependence 
of physical quantities has been formulated. 

Wc have also analyzed the averaged spin-field distribution functions obtained by the high- 
tcmpcrature expansion and the CLS theory for the SK model. Cumulants of the distri- 
bution functions show that the CLS theory gives excellent agreement with simulations at a 
relative low temperature where the high-temperature expansion is not very useful. There- 
fore the CLS theory, if not exact, serves as a reliable tool to analyze dynamics of the SK 
model. 
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Appendix A: Covariances 

We show in this Appendix the procedure to calculate the covariance between two physical quan- 
tities A{a, {Jij}) and B{a, {Jij}) for the SK model defined as 



Cov{AB) = 



AB) 



A) B 



Here ( ) is the thermal average and [ ] represents the configurational average. The exchange 
interaction obeys the Gaussian distribution with vanishing mean and variance J^/N, 

We use the replica method to carry out the sample average 

J2A{cr\{Jij})B{a\{Jij}) 







'{AB) 


= lim 









= lim 







exp[/3^J,,-^<a^) 

\ i<j ce J 

J2Aia\{J.,})B{cT^{J,,}) 
cr 

\ i<j ce ) 



(A-l) 



(A2) 



where replica indices (a = 1, 2, • • • , n) are introduced. 

It is necessary to calculate the average of various powers of exchange interactions to derive the 
expressions of covariances in §4. It is straightforward to show by integration that 

Jabexpl /3Jab^(7^cr^ 

\ Oi J 



^/^^a^^expJ-^'^' 



and 



2N 



Jib exp I 13 Jab < 

V a 



(A-3) 



N 



2N 



(A-4) 



where we have omitted terms which vanish in the n — ^ limit. Higher powers of the interaction 
lead to contributions with higher-power dependence on 1/N. That is, for A; > 3, 

4^feexp( /3J„bE« 

\ a / 
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< 



3J^ 



(A-5) 



Therefore we can ncgicct such higher-order terms in the hmit N ^ oo. 
Let us now work on Cov(i?X)i^i) ^-s an example. Exphcitly written, 



El 



Coy (^H^hi^ = hrn^ 

( E JcdJdef^We ) exp I /3 E -^ij E ^i^l ) 

E(~^E'^«^>^a^M (E-^cd-^deCTcf^e J 
Cr \ ah J \cde J 



(A-6) 



It is straightforward to check that we need only terms in which there are overlaps between two 
interactions. The two quantities H and have four such cases, i.e., ab = cd, ab = dc, ab = de 

and ab = ed. Prom the symmetry under the exchange of site indices, these four cases give identical 
contributions 

Coy ^H^hi^ = 



1 



lim 4 — 



E JabJbe4^l exp /3 E -^ij E ^^j" 



abe 



E JlbJbe<yl'yUl(yl exp ( /3 E Jij E ^i'^i 

abe \ i<j cc / 



The averages over J^^ and J^g are carried out independently. Using eqs. (A-3) and (A-4), the above 
equation is found to be 



lim 



abe \ / \ / cr \ a 



Y.-b->l-Wa-l 1 -P<| ^E( E-r-: 



2N 



The expression in the exponent has been rearranged so that the problem reduces to that of a 
single-site system after the Hubbard-Stratonovitch transformation. The result of spin trace is, in 
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the paramagnetic phase, 



Gov 



= -2NJ^p. (A-7) 
The other covariances can be calculated in the same way. 

Appendix B: Spin-field distribution function 

The dynamical spin-field distribution function is expressed by the characteristic function defined 
in eq. (^.4|). In this Appendix, we give detailed calculations of this function. 

The method of high-temperature series expansion developed in §2 is applicable to the evaluation 
of the characteristic function 



G{k,a) 



1 / ^ 1 -I- fJcr, 



exp ^2 Jij'^j^ 



(B-1) 



The first-order approximation in /3 is obtained in the ordinary way used in equilibrium calculations 
with P replaced by (3eS ■ This procedure is almost the same as the study by Thomsen et aZ.i^ except 
for the factor (1 + aai)/2 



(J 



1 / ^ 1 + (J(Ti 



(B-2) 



1 ' ' ij 

We use the replica method and carry out the average with respect to the distribution of {Jij}, 
so that the function is written as 



G(i)(A:,cj) = lim 



1 1 + craj 

]^EE^r^exp 



<T i 



ik ^ JijCj] \ exp /3cff ^1 '^ji'^ji ) 



1- J ^ -'^ + '^'^J^ I '^^/^cfT 2 , A^ett 1 



ikJ^P, 

w 



ill) 



iim<'lyyi±^ 



exp 



2Ar 



2N 



ikaf 



2 f2 



+ 



'eff 



(B-3) 



(al) \ i 

where = X)o=i '^a^b- Application of the Hubbard-Stratonovitch transformation leads to 



G'-^\k,a) = lim J exp 

n^O 1 



07 > exp 



E'^ 

(07) 
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.7 



> . 



(al) V / ^'^^^ / 

We evaluate the integrals at the saddle point of qa'^: 



(B-4) 



n(n-l)/2 



exp 



(07) > 



E;^ E exp ( ^ ^ 

f i \ (a7)^(al) i 



7 



(al) 



(B-5) 



Furthermore, we assume the system in the paramagnetic phase, that is = ^a,-^. Finally, we 
obtain the first-order characteristic function as 



lim < 



2 t2' 



exp 



n(n-l)/2 



1^1 + aa/ 



^ E^ E '^2""' exp(ifcj2^efr<T,') 



= exp^ 
+ i(7sin(A;J^/3efr)) 



^)i(cOs(fcP/3eff) 



(B-6) 



The higher order approximation needs covarianccs between (1/A^) + exp(i/c^^- JijCj) 

and yl((7, { Jjj}), where ^ stands for the quantities appearing in the expansion of the dynamical 
probability distribution function. The difference from the calculations of other quantities is the 
existence of {ikY^j JijGj) in the exponent. For this reason, the following term is included in 
equations: 



Jib exp I /^eflf^b E + E "^"cO-c 
\ {ij) c 



JL (^^{PesJabUab + ikJab^Tb) 
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exp I l3esJab X] + -^"cO-c I 

V Hi)¥=(ab) c+h ) 



(B.7) 

where the coefficient J^^^^ comes from the quantity A(a^ K-^ij])- The index I is smaller than or equal 
to three in this expansion, because the quantity A contains J^^ at most. 
The following relations are useful for sample averages: 



exp (^f3esJabUab + ikJab^l^ 



exp ( ^(/?efr^a& + ^ikPeScrl'^ab " k'^) 



(B-8) 



Jab exp {j3esJabUab + ikJabC^b) 



J 



X exp ( —{Pesulb + 2ikPeS(^luab - k^) 



(B-9) 



Jab exp (/3eff-'a6^^a6 + ikJabCrfj 



^ iV2 



X exp [ ^(/?eff^^ab + 2ikl3eS(^bUab - k^) 



(B-10) 



In the case of Z = 3, the leading term is of order l/N"^. Therefore there is no contribution in the 
thermodynamic limit and we pay attention only to the case of Z < 2. 

We explain a part of the calculation of the second-order approximation. The second-order term 
is expressed as 



1 + crcTj 



P^h lim { 

n->0 



(T abed 



exp ik X Jdecrl + PeS X] -^y 
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N 



(T abed 



exp ik Y JdeCTe + PcS Y ■^ii'^ij 



(B-ll) 



The overlap between the interactions makes the covariance non-vanishing, but this does not mean 
that all overlaps give finite contributions. In the present case, three configurations d = a,b and c 
have non-vanishing covariances and others are zero. 

Let us consider the configuration d = a as an example. The difference between the first and the 
second terms is only in the replica index of ae in the exponent. Thus we explain only the procedure 
to treat the second term. The procedure is similar to the equilibrium calculation!) and the second 
term is found to be expressed as 



P^bi lim . 



Y JabJacfylcrl exp ( ik ^ JaeCrl + /3eff Y 
CT abc \ e (jj) 



1 



CT abc 



Jab exp I ik V Jabcr'i + fS^sJabU. 



J-ab 



Jac exp \ ik Y JacO-l + PcsJabUa 



exp f ik Y •^aeO'e + PeS Y 

\ e^6,c {ij)^{ab),(ac) 



(B-12) 



J 



^^^^ n^o] iv3^ ^ (^b^^liPcSUab + ikal){PeSUac + ika, 



(J abc 



J' 



X exp — 5^ /?eff4 + E 2^fc/3efrt^, 



2 t2 



k^J 



e^ae 



(B-13) 



Applying the Hubbard-Stratonovitch transformation and the saddle-point method, we obtain in 
the paramagnetic phase 



ATS 



<T abc I 7<5 



-k^alulalal > exp — + ikj'^PeScrl ■ 



(B-14) 



By tracing out the spin configurations, we find the only non- vanishing contribution in the first term 
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for the case 7 = (5 = 1. In this way the second term of eq. (B-11) is written as 



■ cos(A; J^/3efr) 



(B-15) 



in the hmit of ^ 00. 
The first term is obtained by the exchange — > as 



J 



ATS 



<T ahc 75 



\ 7 <5 



-k^alalalal ) exp — + ikj'^PeScrl 



(B-16) 



Non- vanishing contributions come from terms satisfying 7 = 5 = 1: 



f3%{ 



cos (/cJ^ /Jeff) 



-2A: J^/?eff exp — sin(A;P/3eff) 



-A;2j^exp (-'^] sm{kj'M 



(B-17) 



Subtraction of eq. ( B-15 ) from eq. (B-17) gives the covariance of the configuration d = a in eq. (B-12) 
as 



PHA -2A:j4/3effexp ( - ^ ) sm{kJ^pes) 



k^j'^ exp 



k^P 



sin(/cJ^/3efr) 



(B-18) 



For the configurations of d = b and c, the covariances are identical and the result is 

pnA-kJ^pesexp ( 1 sm{kP/3,s)\ . 



(B-19) 



Contributions of all configurations are summed up to yield 



G(2)(A;,a) = P%{ -4A;jVeffexp sm{kJ^P^s 
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2 



sin(fcj2/3eff) >. 



(B-20) 



The third-order term can be calculated similarly. 
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Figure Captions 

Fig. 1. The relaxation of [{H)]/N of the SK model by Monte Carlo simulation, high-temperature 

(3) 

expansions, and the approximation of a Boltzmann form with the inverse temperature . Average 
and standard deviation are shown for 50 samples simulated. The system size is A'^ = 1000 in 
simulations. 

Fig. 2. The relaxation of the expectation value of J2 per spin for the SK model. 

Fig. 3. The relaxation of the expectation value of Y^CFih^ per spin for the SK model. 

Fig. 4. The relaxation of the expectation value of Jijhihj per spin for the SK model. 

Fig. 5. The relaxation of [{H)]/N of the two-dimensional it J model by Monte Carlo simulation, 

(3) 

the first order approximation and the -approximation. Average and standard deviation are 
shown for 50 samples simulated. The system size is A" = 50^ in simulations. 

Fig. 6. The relaxation of the expectation value of X] hf per spin for the two-dimensional zt J 
model. 

Fig. 7. The relaxation of the expectation value of J^CFihf per spin for the two-dimensional it J 
model. 

Fig. 8. The relaxation of the expectation value of J2 Jijhihj per spin for the two-dimensional 
lb J model. 

Fig. 9. Time evolution of the spin-field distribution functions P(/i,ibl) at T = 2. Circles and 
squares are results of simulations. Solid lines denote the theoretical prediction of the third-order 
approximation. 

Fig. 10. Time evolution of physical quantities by the third-order series expansions (full curves) 
and simulations (dotted curves, standard deviations omitted, N = 1000 and 50 samples) when the 
SK model is quenched from Tq = oo to T = 2. 

Fig. 11. Time evolution of correlation coefficients with ^ Jijhihj when the system is quenched 

(3) 

from To = oo to r = 2 obtained by using /^^g . 

Fig. 12. Time evolution of the terms in the dynamical probability distribution function evaluated 
by the third-order series expansions for the SK model quenched from Tq = oo to T = 5. 

Fig. 13. Time evolution of the terms in the dynamical probability distribution function evaluated 
by the third-order series-expansions for the SK model quenched from Tq = oo to T = 2. 

Fig. 14. Time evolution of the terms in the dynamical probability distribution function for the 
SK model quenched from Tq = oo to T = 1 = Tc by simulations. The system size is A^ = 1000. 
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The results are 50-sample averaged. 

Fig. 15. Time evolution of the terms in the dynamical probability distribution function for the 

(3) 

two-dimensional ± J model quenched from Tq = oo to T = 5 by the /S^g -approximations. 

Fig. 16. Time evolution of the terms in the dynamical probability distribution function for the 
two-dimensional ±J model quenched from Tq = oo to T = 2 by simulations. The system size is 
N = 50^ and the results are averaged over 50 samples. 

Fig. 17. First-order cumulants of the spin-field distribution functions P{h, ±1) at T = 2. The 
dotted line is by the high-temperature expansion and the full line is for the CLS theory. Simulation 
results {N = 5000, 100 samples) are also shown for comparison. 

Fig. 18. Second-order cumulants of the spin-field distribution functions P(/i, ±1) at T = 2. 
Symbols are the same as in Fig. 17. 
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